//	--------------------------------------------------------------------
//	Copyright(C) 2006,2007 Zhao Yukun. All rights reserved.
//	This file is a part of the Xeres system.
//	Xeres 2007
//	Contact : xeres.engine@gmail.com
//	$(reservedInfo)
//	$(reservedDeclaration)
//	--------------------------------------------------------------------
///	@file	<matrix.h>
///	@path	~/src/lib/math/
///	@date	2007/09/16
///	@desc	Matrix, 4x4 space matrix represented.

#pragma once

#include "xeres.h"

#include "lib/math/math_ext.h"
#include "lib/math/pair.h"
#include "lib/math/vector.h"

namespace xeres {

	/*!
		\class		Matrix
		\brief		4x4 space matrix.
	*/
	template< typename _Ty > class Matrix : public ArithmeticOp< _Ty >
	{
	public:

		typedef _Ty _Scalar;

		//@ data
		union{

			// linear access
			_Scalar m16[16];

			// easy access
			_Scalar m4x4[4][4];

		};

		// dimension ( accessible member )
		static const size_t DIMENSION = 16;

	public:
		/// \ctor
		Matrix( void )
		{
		}
		/// \ctor_copy
		Matrix( const Matrix& m )
		{
			memcpy( this , &m , sizeof(Matrix) );
		}
		/// \ctor_set
		Matrix(
			const _Scalar t00 , const _Scalar t01 , const _Scalar t02 , const _Scalar t03 ,
			const _Scalar t10 , const _Scalar t11 , const _Scalar t12 , const _Scalar t13 ,
			const _Scalar t20 , const _Scalar t21 , const _Scalar t22 , const _Scalar t23 ,
			const _Scalar t30 , const _Scalar t31 , const _Scalar t32 , const _Scalar t33 )
		{
#	define AS( I , J )	m4x4[I][J] = t##I####J##;
			AS(0,0);AS(0,1);AS(0,2);AS(0,3);
			AS(1,0);AS(1,1);AS(1,2);AS(1,3);
			AS(2,0);AS(2,1);AS(2,2);AS(2,3);
			AS(3,0);AS(3,1);AS(3,2);AS(3,3);
#	undef	AS
		}

		// converter
		template< typename _TOther > Matrix( const Matrix<_TOther>& m )
		{
			for( int i = 0 ; i < 16 ; ++i )
				m16[i] = _Scalar(m.m16[i]);
		}

		// set method ( most not recommended call )
		void set(
			const _Scalar t00 , const _Scalar t01 , const _Scalar t02 , const _Scalar t03 ,
			const _Scalar t10 , const _Scalar t11 , const _Scalar t12 , const _Scalar t13 ,
			const _Scalar t20 , const _Scalar t21 , const _Scalar t22 , const _Scalar t23 ,
			const _Scalar t30 , const _Scalar t31 , const _Scalar t32 , const _Scalar t33 )
		{
#	define AS( I , J )	m4x4[I][J] = t##I####J##;
			AS(0,0);AS(0,1);AS(0,2);AS(0,3);
			AS(1,0);AS(1,1);AS(1,2);AS(1,3);
			AS(2,0);AS(2,1);AS(2,2);AS(2,3);
			AS(3,0);AS(3,1);AS(3,2);AS(3,3);
#	undef	AS
		}

		// set single value
		void single( const _Scalar val )
		{
			m16[ 0] = m16[ 1] = m16[ 2] = m16[ 3] =
			m16[ 4] = m16[ 5] = m16[ 6] = m16[ 7] =
			m16[ 8] = m16[ 9] = m16[10] = m16[11] =
			m16[12] = m16[13] = m16[14] = m16[15] = val;
		}

		// copy from _Tying point array
		void copy( const _Scalar * f )
		{
			memcpy( m16 , f , sizeof(_Scalar) * DIMENSION );
		}

		// set zero matrix
		void set_zero( void )
		{
			memset( m16 , 0 , sizeof(_Scalar) * DIMENSION );
		}

		// load identity matrix
		void identity( void )
		{
			static const _Scalar ALIGN16 idmm[] = {
				1,0,0,0,
				0,1,0,0,
				0,0,1,0,
				0,0,0,1,
			};
			memcpy( this , idmm , sizeof(_Scalar) * DIMENSION );
		};

		// get static identity
		static const Matrix<_Scalar>& _identity( void )
		{
			static const ALIGN16 Matrix<_Scalar> idmm(
				1,0,0,0,
				0,1,0,0,
				0,0,1,0,
				0,0,0,1
			);
			return idmm;
		}

	public:

		// self transpose
		void transpose( void )
		{
			register _Scalar t;
#	define swap( a , b ) { t = a; a = b; b = t; }
			swap( m4x4[1][0] , m4x4[0][1] );
			swap( m4x4[2][0] , m4x4[0][2] );
			swap( m4x4[2][1] , m4x4[1][2] );
			swap( m4x4[3][0] , m4x4[0][3] );
			swap( m4x4[3][1] , m4x4[1][3] );
			swap( m4x4[3][2] , m4x4[2][3] );
#	undef	swap
		}

		// make transpose
		Matrix transposed( void ) const
		{
			Matrix out;
			out.m4x4[0][0]	= m4x4[0][0];
			out.m4x4[1][0]	= m4x4[0][1];
			out.m4x4[2][0]	= m4x4[0][2];
			out.m4x4[3][0]	= m4x4[0][3];
			out.m4x4[0][1]	= m4x4[1][0];
			out.m4x4[1][1]	= m4x4[1][1];
			out.m4x4[2][1]	= m4x4[1][2];
			out.m4x4[3][1]	= m4x4[1][3];
			out.m4x4[0][2]	= m4x4[2][0];
			out.m4x4[1][2]	= m4x4[2][1];
			out.m4x4[2][2]	= m4x4[2][2];
			out.m4x4[3][2]	= m4x4[2][3];
			out.m4x4[0][3]	= m4x4[3][0];
			out.m4x4[1][3]	= m4x4[3][1];
			out.m4x4[2][3]	= m4x4[3][2];
			out.m4x4[3][3]	= m4x4[3][3];
			return out;
		}

		// self inverse
		void inverse( void )
		{
			Matrix<_Scalar> tmp;
			_Scalar det = _Scalar( det4x4( m16 ) );
			if( ! xeres::is_zero( det ) )
			{
				adjointx<_Scalar>( tmp.m16 , m16 );
				if( !is_equal<_Scalar>( det , 1 ) )
				{
					det = recip( det );
					for( int i = 0 ; i < 16 ; ++i )
						tmp.m16[i] *= det;
					*this = tmp;
				}
			}
			else
			{
				identity();
			}
		}

		// make invert
		Matrix invert( void ) const
		{
			Matrix out;
			_Scalar det = _Scalar( det4x4( m16 ) );
			if( ! xeres::is_zero( det ) )
			{
				adjointx<_Scalar>( out.m16 , m16 );
				if( !is_equal<_Scalar>( det , 1 ) )
				{
					det = recip( det );
					for( int i = 0 ; i < 16 ; ++i )
						out.m16[i] *= det;
				}
			}
			else
			{
				out.identity();
			}
			return out;
		}

		// be the base of a coordinate system
		void be_base( const Vector<_Scalar>& e1 , const Vector<_Scalar>& e2 , const Vector<_Scalar>& e3 )
		{
			static const _Scalar endline[] = { 0 , 0 , 0 , 0 , 1 };
			//
			m16[0 ] = e1[0];
			m16[1 ] = e1[1];
			m16[2 ] = e1[2];
			m16[3 ] = 0;
			//
			m16[4 ] = e2[0];
			m16[5 ] = e2[1];
			m16[6 ] = e2[2];
			m16[7 ] = 0;
			//
			m16[8 ] = e3[0];
			m16[9 ] = e3[1];
			m16[10] = e3[2];
			//
			memcpy( &m16[11] , endline , sizeof(endline) );
		}

        // set a translation matrix
		void be_translation( const Vector<_Scalar>& t )
		{
			identity();
			m16[12] = t[0];
			m16[13] = t[1];
			m16[14] = t[2];
		}

		void be_translation( const _Scalar tx , const _Scalar ty , const _Scalar tz )
		{
			identity();
			m16[12] = tx;
			m16[13] = ty;
			m16[14] = tz;
		}

        // set a rotation matrix
		void be_rotation( const _Scalar ax , const _Scalar ay , const _Scalar az , const _Scalar n_radians )
		{
			Vector<_Scalar> axis( ax , ay , az );
			axis.normalize();

			_Scalar sine , cosine;
			sincos( n_radians , &sine , &cosine );

			_Scalar ab = axis[0] * axis[1] * (1 - cosine);
			_Scalar bc = axis[1] * axis[2] * (1 - cosine);
			_Scalar ca = axis[2] * axis[0] * (1 - cosine);

			_Scalar t;

			identity();

			t = axis[0] * axis[0];
			m4x4[0][0] = t + cosine * (1 - t);
			m4x4[2][1] = bc - axis[0] * sine;
			m4x4[1][2] = bc + axis[0] * sine;

			t = axis[1] * axis[1];
			m4x4[1][1] = t + cosine * (1 - t);
			m4x4[2][0] = ca + axis[1] * sine;
			m4x4[0][2] = ca - axis[1] * sine;

			t = axis[2] * axis[2];
			m4x4[2][2] = t + cosine * (1 - t);
			m4x4[1][0] = ab - axis[2] * sine;
			m4x4[0][1] = ab + axis[2] * sine;
		}

		void be_rotation( const Vector<_Scalar>& ax , const _Scalar n_radians )
		{
			Vector<_Scalar> axis = ax.normalized();

			_Scalar sine , cosine;
			sincos( n_radians , &sine , &cosine );

			_Scalar ab = axis[0] * axis[1] * (1 - cosine);
			_Scalar bc = axis[1] * axis[2] * (1 - cosine);
			_Scalar ca = axis[2] * axis[0] * (1 - cosine);

			float t;

			identity();

			t = axis[0] * axis[0];
			m4x4[0][0] = t + cosine * (1 - t);
			m4x4[2][1] = bc - axis[0] * sine;
			m4x4[1][2] = bc + axis[0] * sine;

			t = axis[1] * axis[1];
			m4x4[1][1] = t + cosine * (1 - t);
			m4x4[2][0] = ca + axis[1] * sine;
			m4x4[0][2] = ca - axis[1] * sine;

			t = axis[2] * axis[2];
			m4x4[2][2] = t + cosine * (1 - t);
			m4x4[1][0] = ab - axis[2] * sine;
			m4x4[0][1] = ab + axis[2] * sine;
		}

		// set a scale matrix
		void be_scale( const _Scalar sx , const _Scalar sy , const _Scalar sz )
		{
			identity();
			m16[0 ] = sx;
			m16[5 ] = sy;
			m16[10] = sz;
		}
		void be_scale( const Vector<_Scalar>& s )
		{
			identity();
			m16[0 ] = s[0];
			m16[5 ] = s[1];
			m16[10] = s[2];
		}

		// multiply 2 matrices
		template< typename _T1 , typename _T2 , typename _T3 > static inline void mat_multiply(
			const Matrix<_T1>& a , const Matrix<_T2>& b , Matrix<_T3>& c )
		{
			typedef double _TT;

			register const _T1 * ap = &a.m4x4[0][0];
			register const _T2 * bp = &b.m4x4[0][0];
			register       _T3 * cp = &c.m4x4[0][0];

			register _TT a0, a1, a2, a3;

			a0 = _TT(ap[0]);
			a1 = _TT(ap[1]);
			a2 = _TT(ap[2]);
			a3 = _TT(ap[3]);

			cp[0] = _Ty( a0 * bp[0]  + a1 * bp[4]  + a2 * bp[8 ] + a3 * bp[12] );
			cp[1] = _Ty( a0 * bp[1]  + a1 * bp[5]  + a2 * bp[9 ] + a3 * bp[13] );
			cp[2] = _Ty( a0 * bp[2]  + a1 * bp[6]  + a2 * bp[10] + a3 * bp[14] );
			cp[3] = _Ty( a0 * bp[3]  + a1 * bp[7]  + a2 * bp[11] + a3 * bp[15] );

			a0 = _TT(ap[4]);
			a1 = _TT(ap[5]);
			a2 = _TT(ap[6]);
			a3 = _TT(ap[7]);

			cp[4] = _Ty( a0 * bp[0]  + a1 * bp[4]  + a2 * bp[8 ] + a3 * bp[12] );
			cp[5] = _Ty( a0 * bp[1]  + a1 * bp[5]  + a2 * bp[9 ] + a3 * bp[13] );
			cp[6] = _Ty( a0 * bp[2]  + a1 * bp[6]  + a2 * bp[10] + a3 * bp[14] );
			cp[7] = _Ty( a0 * bp[3]  + a1 * bp[7]  + a2 * bp[11] + a3 * bp[15] );

			a0 = _TT(ap[8 ]);
			a1 = _TT(ap[9 ]);
			a2 = _TT(ap[10]);
			a3 = _TT(ap[11]);

			cp[8 ] = _Ty( a0 * bp[0]  + a1 * bp[4]  + a2 * bp[8 ] + a3 * bp[12] );
			cp[9 ] = _Ty( a0 * bp[1]  + a1 * bp[5]  + a2 * bp[9 ] + a3 * bp[13] );
			cp[10] = _Ty( a0 * bp[2]  + a1 * bp[6]  + a2 * bp[10] + a3 * bp[14] );
			cp[11] = _Ty( a0 * bp[3]  + a1 * bp[7]  + a2 * bp[11] + a3 * bp[15] );

			a0 = _TT(ap[12]);
			a1 = _TT(ap[13]);
			a2 = _TT(ap[14]);
			a3 = _TT(ap[15]);

			cp[12] = _Ty( a0 * bp[0]  + a1 * bp[4]  + a2 * bp[8 ] + a3 * bp[12] );
			cp[13] = _Ty( a0 * bp[1]  + a1 * bp[5]  + a2 * bp[9 ] + a3 * bp[13] );
			cp[14] = _Ty( a0 * bp[2]  + a1 * bp[6]  + a2 * bp[10] + a3 * bp[14] );
			cp[15] = _Ty( a0 * bp[3]  + a1 * bp[7]  + a2 * bp[11] + a3 * bp[15] );
		}

		// multiply 2 matrix
		static inline Matrix mat_multiply( const Matrix& a , const Matrix& b )
		{
			Matrix tmp;
			mat_multiply( a , b , tmp );
			return tmp;
		}
		// this = this * m
		void multiply_to( const Matrix& m )
		{
			Matrix tmp;
			mat_multiply( *this , m , tmp );
			*this = tmp;
		}
		// this = m * this
		void multiply_from( const Matrix& m )
		{
			Matrix tmp;
			mat_multiply( m , *this , tmp );
			*this = tmp;
		}

		// translate matrix
		void translate( const Vector<_Scalar>& t )
		{
			m4x4[3][0] += t.x * m4x4[0][0] + t.y * m4x4[1][0] + t.z * m4x4[2][0];
			m4x4[3][1] += t.x * m4x4[0][1] + t.y * m4x4[1][1] + t.z * m4x4[2][1];
			m4x4[3][2] += t.x * m4x4[0][2] + t.y * m4x4[1][2] + t.z * m4x4[2][2];
			m4x4[3][3] += t.x * m4x4[0][3] + t.y * m4x4[1][3] + t.z * m4x4[2][3];
		}
		void translate( const _Scalar tx , const _Scalar ty , const _Scalar tz )
		{
			m4x4[3][0] += tx * m4x4[0][0] + ty * m4x4[1][0] + tz * m4x4[2][0];
			m4x4[3][1] += tx * m4x4[0][1] + ty * m4x4[1][1] + tz * m4x4[2][1];
			m4x4[3][2] += tx * m4x4[0][2] + ty * m4x4[1][2] + tz * m4x4[2][2];
			m4x4[3][3] += tx * m4x4[0][3] + ty * m4x4[1][3] + tz * m4x4[2][3];
		}

		// rotate matrix
		void rotate( const _Scalar ax , const _Scalar ay , const _Scalar az , const _Scalar n_radians )
		{
			Matrix m;
			m.be_rotation( ax , ay , az , n_radians );
			multiply_to(m);
		}
		void rotate( const Vector<_Scalar>& axis , const _Scalar n_radians )
		{
			Matrix m;
			m.be_rotation( axis , n_radians );
			multiply_to(m);
		}

		// scale matrix
		void scale( const _Scalar sx , const _Scalar sy , const _Scalar sz )
		{
			m4x4[0][0] = (sx * m4x4[0][0]);
			m4x4[0][1] = (sx * m4x4[0][1]);
			m4x4[0][2] = (sx * m4x4[0][2]);
			m4x4[0][3] = (sx * m4x4[0][3]);

			m4x4[1][0] = (sy * m4x4[1][0]);
			m4x4[1][1] = (sy * m4x4[1][1]);
			m4x4[1][2] = (sy * m4x4[1][2]);
			m4x4[1][3] = (sy * m4x4[1][3]);

			m4x4[2][0] = (sz * m4x4[2][0]);
			m4x4[2][1] = (sz * m4x4[2][1]);
			m4x4[2][2] = (sz * m4x4[2][2]);
			m4x4[2][3] = (sz * m4x4[2][3]);
		}
		void scale( const Vector<_Scalar>& s )
		{
			m4x4[0][0] = (s[0] * m4x4[0][0]);
			m4x4[0][1] = (s[0] * m4x4[0][1]);
			m4x4[0][2] = (s[0] * m4x4[0][2]);
			m4x4[0][3] = (s[0] * m4x4[0][3]);

			m4x4[1][0] = (s[1] * m4x4[1][0]);
			m4x4[1][1] = (s[1] * m4x4[1][1]);
			m4x4[1][2] = (s[1] * m4x4[1][2]);
			m4x4[1][3] = (s[1] * m4x4[1][3]);

			m4x4[2][0] = (s[2] * m4x4[2][0]);
			m4x4[2][1] = (s[2] * m4x4[2][1]);
			m4x4[2][2] = (s[2] * m4x4[2][2]);
			m4x4[2][3] = (s[2] * m4x4[2][3]);
		}

		// origin
		Vector<float> origin( void ) const
		{
			return *reinterpret_cast<const Vector<float>*>( m4x4[3] );
		}

		// axis
		const Vector<float>& axis( int a ) const
		{
			XS_ASSERT( a >= 0 && a < 3 );
			return *reinterpret_cast<const Vector<float>*>( m4x4[a] );
		}

		// xform methods

		template<typename _T1 , typename _T2> void xformp( const Vector<_T1>& from , Vector<_T2>& to ) const
		{
			_Scalar x = m4x4[0][0] * from[0] + m4x4[1][0] * from[1] + m4x4[2][0] * from[2] + m4x4[3][0];
			_Scalar y = m4x4[0][1] * from[0] + m4x4[1][1] * from[1] + m4x4[2][1] * from[2] + m4x4[3][1];
			_Scalar z = m4x4[0][2] * from[0] + m4x4[1][2] * from[1] + m4x4[2][2] * from[2] + m4x4[3][2];
			_Scalar w = m4x4[0][3] * from[0] + m4x4[1][3] * from[1] + m4x4[2][3] * from[2] + m4x4[3][3];
			if( !is_equal( w , _Scalar(1) ) )
			{
				w = recip( w );
				x *= w;
				y *= w;
				z *= w;
			}
			to.set( x , y , z );
		}
		template<typename _T1 , typename _T2> void xformv( const Vector<_T1>& from , Vector<_T2>& to ) const
		{
			_Scalar x = m4x4[0][0] * from[0] + m4x4[1][0] * from[1] + m4x4[2][0] * from[2];
			_Scalar y = m4x4[0][1] * from[0] + m4x4[1][1] * from[1] + m4x4[2][1] * from[2];
			_Scalar z = m4x4[0][2] * from[0] + m4x4[1][2] * from[1] + m4x4[2][2] * from[2];
			to.set( x , y , z );
		}
		template<typename _T1 , typename _T2> void xformn( const Vector<_T1>& from , Vector<_T2>& to ) const
		{
			_Scalar x = m4x4[0][0] * from[0] + m4x4[0][1] * from[1] + m4x4[0][2] * from[2];
			_Scalar y = m4x4[1][0] * from[0] + m4x4[1][1] * from[1] + m4x4[1][2] * from[2];
			_Scalar z = m4x4[2][0] * from[0] + m4x4[2][1] * from[1] + m4x4[2][2] * from[2];
			to.set( x , y , z );
		}

		// xform in flip order

		template<typename _T1 , typename _T2> void xformp_flip( const Vector<_T1>& from , Vector<_T2>& to ) const
		{
			_Scalar x = m4x4[0][0] * from[0] + m4x4[0][1] * from[1] + m4x4[0][2] * from[2] + m4x4[0][3];
			_Scalar y = m4x4[1][0] * from[0] + m4x4[1][1] * from[1] + m4x4[1][2] * from[2] + m4x4[1][3];
			_Scalar z = m4x4[1][0] * from[0] + m4x4[2][1] * from[1] + m4x4[2][2] * from[2] + m4x4[2][3];
			_Scalar w = m4x4[3][0] * from[0] + m4x4[3][1] * from[1] + m4x4[3][2] * from[2] + m4x4[3][3];
			if( !is_equal( w , _Scalar(1) ) )
			{
				w = recip( w );
				x *= w;
				y *= w;
				z *= w;
			}
			to.set( x , y , z );
		}
		template<typename _T1 , typename _T2> void xformv_flip( const Vector<_T1>& from , Vector<_T2>& to ) const
		{
			_Scalar x = m4x4[0][0] * from[0] + m4x4[0][1] * from[1] + m4x4[0][2] * from[2];
			_Scalar y = m4x4[1][0] * from[0] + m4x4[1][1] * from[1] + m4x4[1][2] * from[2];
			_Scalar z = m4x4[2][0] * from[0] + m4x4[2][1] * from[1] + m4x4[2][2] * from[2];
			to.set( x , y , z );
		}
		template<typename _T1 , typename _T2> void xformn_flip( const Vector<_T1>& from , Vector<_T2>& to ) const
		{
			_Scalar x = m4x4[0][0] * from[0] + m4x4[1][0] * from[1] + m4x4[2][0] * from[2];
			_Scalar y = m4x4[0][1] * from[0] + m4x4[1][1] * from[1] + m4x4[2][1] * from[2];
			_Scalar z = m4x4[0][2] * from[0] + m4x4[1][2] * from[1] + m4x4[2][2] * from[2];
			to.set( x , y , z );
		}

		// comparison
		bool is_identity( void ) const
		{
			static const _Scalar ALIGN16 idmm[] = {
				1,0,0,0,
				0,1,0,0,
				0,0,1,0,
				0,0,0,1, };
#	define CMP(X)	if( !is_equal( m16[X] , idmm[X] ) ) return false;
			CMP(0);CMP(1);CMP(2);CMP(3);
			CMP(4);CMP(5);CMP(6);CMP(7);
			CMP(8);CMP(9);CMP(10);CMP(11);
			CMP(12);CMP(13);CMP(14);CMP(15);
			return true;
#	undef	CMP
		}
		bool is_zero( void ) const
		{
#	define CMP(X)	if( !is_zero( m16[X] ) ) return false;
			CMP(0);CMP(1);CMP(2);CMP(3);
			CMP(4);CMP(5);CMP(6);CMP(7);
			CMP(8);CMP(9);CMP(10);CMP(11);
			CMP(12);CMP(13);CMP(14);CMP(15);
			return true;
#	undef	CMP
		}

		bool operator == ( const Matrix& m ) const
		{
#	define CMP(X)	if( !is_equal( m16[X] , m.m16[X] ) ) return false;
			CMP(0);CMP(1);CMP(2);CMP(3);
			CMP(4);CMP(5);CMP(6);CMP(7);
			CMP(8);CMP(9);CMP(10);CMP(11);
			CMP(12);CMP(13);CMP(14);CMP(15);
			return true;
#	undef	CMP
		}
		bool operator != ( const Matrix& m ) const
		{
#	define CMP(X)	if( !is_equal( m16[X] , m.m16[X] ) ) return true;
			CMP(0);CMP(1);CMP(2);CMP(3);
			CMP(4);CMP(5);CMP(6);CMP(7);
			CMP(8);CMP(9);CMP(10);CMP(11);
			CMP(12);CMP(13);CMP(14);CMP(15);
			return false;
#	undef	CMP
		}

		// index
		const _Scalar& index( int i ) const
		{
			assert( ( i >= 0 ) && ( i < 16 ) );
			return m16[i];
		}

		// index
		_Scalar& index( int i )
		{
			assert( ( i >= 0 ) && ( i < 16 ) );
			return m16[i];
		}

		// element
		const _Scalar& element( int row , int col ) const
		{
			assert( row >= 0 && row < 4 && col >= 0 && col < 4 );
			return m4x4[ row ][ col ];
		}

		// element
		_Scalar& element( int row , int col )
		{
			assert( row >= 0 && row < 4 && col >= 0 && col < 4 );
			return m4x4[ row ][ col ];
		}

		// element
		const _Scalar * operator [] ( int row ) const
		{
			assert( row >= 0 && row < 4 );
			return m4x4[ row ];
		}

		// element
		_Scalar * operator [] ( int row )
		{
			assert( row >= 0 && row < 4 );
			return m4x4[ row ];
		}

		// determinant
		_Scalar det( void ) const
		{
			return _Scalar(det4x4(m16));
		}

		Matrix& operator = ( const Matrix& m )
		{
			memcpy( this , m.m16 , sizeof(_Scalar) * DIMENSION );
			return *this;
		}
		template< typename _Other > const Matrix& operator = ( const Matrix<_Other>& m )
		{
			for( size_t i = 0 ; i < 16 ; ++i )
				m16[i] = _Scalar(m.m16[i]);
			return *this;
		}
	};

} // namespace xeres
